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, Invariant parameterization schemes for the eddy-vorticity flux in the barotropic vorticity equation 

on the beta-plane are constructed and then apphed to turbulence modeling. This construction 
, is realized by the exhaustive description of differential invariants for the maximal Lie invariance 

pseudogroup of this equation using the method of moving frames, which includes finding functional 
bases of differential invariants of arbitrary order, a minimal generating set of differential invari- 
ants and a basis of operators of invariant differentiation in an explicit form. Special attention is 
paid to the problem of two-dimensional turbulence on the beta-plane. It is shown that classical 
hyperdiffusion as used to initiate the energy-enstrophy cascades violates the symmetries of the 
vorticity equation. Invariant but nonlinear hyperdiffusion-like terms of new types are introduced 
and then used in the course of numerically integrating the vorticity equation and carrying out 
pi ^ freely decaying turbulence tests. It is found that the invariant hyperdiffusion scheme is close to 

but not exactly reproducing the k~^ shape of energy spectrum in the enstrophy inertial range. By 
presenting conservative invariant hyperdiffusion terms, we also demonstrate that the concepts of 
■ invariant and conservative parameterizations are consistent. 



^ . 1 Introduction 
> 

As atmospheric and oceanic numerical models get increasingly complex, it becomes more and 
On I more challenging to propose valuable conceptual paradigms for those processes that the model 

is still not able to capture owing to its limited spatial and temporal resolution. This problem is 
, common to all numerical models irrespectively of their eventual degree of sophistication [50, 51]. 

In the beginning of numerical modeling in geophysical fluid dynamics, it was often the lack 
of computer power that dictated which processes had to be parameterized, even in case one 
might already had a rather concise understanding of these processes. As computers became 
more capable, the problem of parameterization shifted to processes occurring on rather fine 
scales where it can be difficult to retrieve accurate experimental data. Accordingly, for various 
5^ . processes that should be taken into account in order to improve the forecast range of a numerical 

model, there is still no satisfactory general understanding. This naturally makes it difficult to 
set up valuable parameterization schemes, which for this reason is usually an elaborate task. 

One the other hand, processes that occur in geophysical fluid dynamics and that can be de- 
scribed using differential equations also might have certain structural or geometrical properties. 
Such properties can be conservation of mass or energy or other fundamental conservation laws. 
Real-world processes are generally also invariant under specific transformation groups, as e.g. 
the Galilean group. This is why one can ask the question whether it is reasonable to construct 
parameterization schemes for processes possessing certain structural features in a manner such 
that these features are preserved in the closed differential equation. In this way, even if a model 
is not able to explicitly resolve processes, loosely speaking, it takes into account some of their 
significant properties. This study was initiated in [32] for the problem of finding invariant tur- 
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bulence closure schemes for the filtered Navier-Stokes equations. In the present paper we aim to 
give a further instance for invariant parameterization schemes by constructing closure ansatzes 
that retain certain Lie point symmetries of the barotropic vorticity equation on the beta-plane. 

This possible stream of constructing geometrically motivated parameterization schemes in 
some sense parallels the present general trend in numerical modeling to design specially adapted 
discretizations of differential equations that capture a range of their qualitative or global features, 
such as conservation laws, a Hamiltonian form or symmetry properties. Especially the possibility 
of constructing discretization schemes that have the same symmetries as the original differential 
equations they are a model of, as proposed and discussed in e.g. [17, 25, 26, 44, 52], is of 
immediate relevance to the present work. This is because, strictly speaking, a discretization 
of a system of differential equation is in practice not enough to set up a valuable numerical 
model. There always has to be a model for the unresolved parts of the dynamics. (Neglecting 
them is in general not a good idea as for nonlinear differential equations these parts will, sooner 
or later, spoil the numerical integration.) Then, if one aims to construct an invariant discrete 
counterpart of some relevant physical model, care should also be taken about the invariance 
characteristics of the processes that are not explicitly resolved. This is where the method 
of finding invariant parameterization schemes comes into play. The combination of invariant 
discretization schemes for the resolved part of the model with invariant parameterization schemes 
for the unresolved parts yields a completely invariant numerical description of a given system of 
differential equations. Such a fully invariant model might be closer to a true geometric numerical 
integration scheme than solely a symmetry preserving discretization without any closure for the 
subgrid-scale terms or with some noninvariant closure. 

Perhaps the most relevant usage of the barotropic vorticity equation is related to two- 
dimensional turbulence. Turbulence on the beta-plane (or, more general, on the rotating sphere) 
is peculiar in that it allows for the combination of turbulent and wave-like effects. It is believed to 
explain the emergence of strong jets and band-like structure on giant planets in our solar system 
and is therefore the subject of intensive investigation, see e.g. [23, 31, 45, 47, 53] and references 
therein. In the present paper we focus on freely decaying turbulence on the beta-plane by using 
invariant hyperdiffusion terms to initiate the energy-enstrophy cascades. These cascades are 
most probable responsible for to the emergence of coherent, stable structures (vortices) that are 
ubiquitous in large-scale geophysical fluid dynamics. 

The outline of the paper is as follows. In the subsequent Section 2, we discuss and slightly 
extend the concept of invariant parameterization schemes as introduced in [32] and [41]. Special 
attention will be paid to methods related to invariant parameterization schemes and inverse 
group classification. In Section 3 we present the maximal Lie invariance algebra and the maximal 
Lie invariance pseudogroup of the barotropic vorticity equation on the beta-plane. A concise 
description of the general method for computing differential invariants of Lie (pseudo)groups 
using the method of moving frames is given in Section 4. In Section 5 the algebra of differential 
invariants is determined for the maximal Lie invariance pseudogroup of the vorticity equation. 
Two examples for invariant parameterization schemes constructed out of existing schemes using 
the invariantization process are presented in Section 6. Section 7 is devoted to the application 
of differential invariants in turbulence on the beta-plane. In particular, invariant hyperdiffusion 
schemes are introduced. The vorticity equation on the beta-plane is integrated numerically 
using both invariant and noninvariant hyperdiffusion and the corresponding energy spectra are 
obtained. In Section 8 we discuss the possibility of deriving invariant parameterization schemes 
that also respect conservation laws. As an example, an invariant diffusion term is constructed 
that preserves the entire maximal Lie invariance pseudogroup of the vorticity equation and 
also preserves conservation of energy, circulation and momentum. The results are summarized 
and further discussed in the final Section 9, in which we also indicate possible future research 
directions in the field of invariant parameterization theory. 
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2 Invariant parameterization schemes 



The problem of finding parameterization or closure schemes for subgrid-scale terms in averaged 
differential equations that admit Lie symmetries of the original (unaveraged) differential equation 
was first raised in [32], see also [33, 43]. Recently, we put this idea into the framework of 
group classification [41] , by showing that any problem of constructing invariant parameterization 
schemes amounts in solving a (possibly complicated) group classification problem. 

As for the classical group classification, there are two principal ways to construct parame- 
terization schemes, the direct and the inverse one [41]. In the direct approach, one replaces 
the terms to be parameterized with arbitrary functions depending on the mean variables and 
derivatives thereof. This is in the line with the general definition of all physical parameterization 
schemes, which are concerned to express the unknown subgrid-scale terms using the information 
included in the grid-scale (mean) quantities. The form of dependency of the arbitrary func- 
tions on the mean variables is guided by physical intuition. It determines the properties of all 
the families of invariant parameterization schemes that can be derived (e.g. the highest order 
of derivatives that can arise). Once the general form of the arbitrary function is chosen, one 
is left with a possibly rather general class of differential equations, which is amendable with 
tools from usual group classification, see e.g. [8, 18, 42]. This in particular will lead to a list of 
families of mutually inequivalent parameterization schemes that admit different Lie invariance 
algebras. One can then select those families that preserve the most essential symmetry features 
of the process to be parameterized. The final step is to suitably narrow the selected families by 
including other desirable physical properties into the invariant parameterization scheme. 

In the present paper, however, we will be solely concerned with the inverse approach, which 
is why we will discuss it in a more extended manner. The inverse approach rests on the fact 
that any system of differential equations can be rewritten in terms of differential invariants of its 
maximal Lie invariance group, provided that the prolongation of the group to the corresponding 
jet space acts semi-regularly [35]. This property can be used in the course of the parameterization 
problem in the following way: Suppose that we are given a Lie group G regarded as important 
to be preserved for valuable parameterization schemes as a Lie symmetry group. Computing 
a basis of differential invariants of G along with a complete set of its independent operators 
of invariant differentiation, see e.g. [19, 20, 36, 39], serves to exhaustively describe the entire 
algebra of differential invariants of G. As any combination of these differential invariants will 
necessarily be invariant with respect to G, assembling them together to a parameterization will 
immediately lead to a closure scheme admitting G as a Lie symmetry group. 

The key question hence lies in the correct selection of a suitable symmetry group. The initial 
point for the selection is given by symmetry properties of the model to be parameterized. In the 
course of the parameterization one can intend to preserve the whole Lie symmetry group of the 
initial model or its proper subgroups. The choice for an invariance group for parameterization 
obviously should not solely be justified using mathematical arguments. Sometimes, it can be 
motivated from obvious physical reasons. If the process to be parameterized can be described 
within the framework of classical mechanics then any reasonable parameterization for that pro- 
cess should be invariant under the Galilean group. Moreover, for turbulence closure schemes, 
scale invariance might be of particular importance. For processes that can be described in the 
framework of a variational principle and respect certain conservation laws, it might be reasonable 
that the parameterization scheme to be developed respects the associated Noether symmetries. 

There are several processes in fiuid mechanics that are intimately linked to the presence of 
certain boundary conditions (e.g., turbulence near walls, boundary layer convection, etc.). For 
such processes the inclusion of the particular boundaries is an integral part in the formulation of 
a parameterization scheme. At first glance, to find invariant parameterization schemes for such 
processes it is inevitable to single out those subgroups of the maximal Lie invariance group G 
of the system C of differential equations describing the process of interest that are compatible 
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with a particular boundary value problem. The main complication with this approach is that 
most of boundary value problems admit no symmetries, see e.g. [12]. At the same time, it is 
more natural to assume that symmetries of C act as equivalence transformations on a joint 
class of physically relevant boundary value problems for £, i.e., these transformations send a 
particular boundary value problem to another problem from the same class. Even the basic 
physical symmetries including shifts in space and time, rotations, scalings, Galilean boosts or 
Lorentz transformations, which are related to fundamental properties of the space and the time 
(homogeneity, isotropy, similarity, Galilean or special relativity principle, respectively), usually 
act on boundary value problems in much the same way as equivalence transformations. This 
is why it is the generation of a group of well-defined equivalence transformations on a properly 
chosen class of boundary value problems that can serve as a criterion for selecting a subgroup 
of G to be taken into account in the course of invariant parameterization of C. 

Employing techniques of inverse group classification does not automatically lead to ready-to- 
use parameterizations, but it gives a frame in which parameterizations can be defined without 
the violation of basic invariance properties. Examples of the violation have been reported in 
the recent literature. See, e.g., [32] for a discussion about the Smagorinsky model in the filtered 
Navier-Stokes equations violating scale invariance and [41] for a note on the Kuo convection 
schemes that describes a Galilean invariant process in a noninvariant fashion. The construction 
of parameterization schemes that fail to preserve essential symmetries can be easily avoided 
by applying the above methods of inverse group classification. This may help to restrict the 
large number of possible closure schemes using geometrical reasoning and thereby may assist in 
finding a proper description of unresolved subgrid-scale processes. 

There is yet a second possibility to construct invariant parameterization schemes using the 
inverse group classification approach, which has not been reported in [41]. It rests on the 
construction of moving frames for the Lie group G with respect to which parameterizations under 
study should be invariant. It is a general feature of a moving frame that it allows constructing 
of invariant counterparts of differential functions. This property enables the construction of an 
invariant parameterization scheme out of a noninvariant one. It is simply necessary to apply 
the moving frame corresponding to the selected Lie group G to the specific closed differential 
equation. More precisely, consider a system C of differential equations 

L\x,U(n)) = ^, l = l,...,m. 

The dependent variables u can be represented according to u = u + u' , where u is the averaged or 
filtered part of the dynamics (i.e. the resolved or grid-scale part) and u' denotes the departure of 
u from the mean or filtered part u (i.e. the subgrid-scale part). Numerical models in geophysical 
fluid dynamics are formulated as equations for the resolved part, which are obtained from = 
by averaging or filtering, leading to 

L'-{x,ui^n),w) =0, / = l,...,m, 

where U are smooth differential functions of their arguments. The particular form of U depends 
on the actual averaging rule chosen and the form of the initial system C. The unknown subgrid- 
scale terms that arise in the course of averaging (e.g. by using the Reynolds averaging rule for 
products, ab = ab + a'b') are collected in the tuple w. These terms have to be parameterized 
in order to close the above system of averaged differential equations. A local parameterization 
scheme assumes a particular functional relation 

W = 9{x, U(^r)) 

between the subgrid-scale and grid-scale quantities. Let there be given a moving frame p'--'^ of 
order j = max(r, n) for the selected Lie group G, see Section 4. Any particular parameterization 
scheme can then be invariantized via replacing U and 9 by their invariantized counterparts, 

t{U{x, ■U(„), w)) = U{p^^^ ■ X, p^^^ ■ 'U(„), t«) and l{9{x, U(r))) = 9{p^^^ ■ x, p^^^ ■ U(^r))- 
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3 Lie symmetries of the vorticity equation 



The barotropic vorticity equation on the beta-plane is a simple but still genuine meteorological 
model. It has the form 

Ct + tpxCy - ipyCx + f^ipx = 0, or r]t + Jiip, ri) = 0. (1) 

Here J {a, b) := Uxby — Uybx, ijj = ip{t, x, y) is the stream function, Q = ipxx + fpyy is the vorticity 
and 7y = C + f = C + fo + /3yis the absolute vorticity under the /3-plane approximation f = fo + /3y 
of the Coriolis parameter f, /? is a nonzero constant parameter (the differential rotation). The 
constant fo is dynamically inessential and can be neglected. 

The maximal Lie invariance algebra 31 of Eq. (1) is spanned by the vector fields 

T> = tdt-xdx-ydy-3ipd^, dt, dy, X{f) = f{t)dx - ft{t)yd^ Z{g) = g{t)d^, 

where the parameters / and g run through the set of smooth functions of i [9, 24]. The vorticity 
equation (1) also admits two independent discrete symmetries, which alternate signs of the pairs 
{t,x) and {y,ip), see [10] for more details. Such discrete symmetries will not be taken into 
account in the course of construction of differential invariants. Any nonzero value of /3 can be 
gauged to one by a scaling transformation. 

The one-parameter Lie (pseudo) groups generated by the above vector fields read 





{t,x,y,ip) {e^'^t,e ^'^x,e 


y,e 




{t,x,y,^p) ^ {t + e2,x,y,ip) 






{t,x,y,^p) ^ {t,x,y + e3,'tp) 






{t, X, y, i)) ^ {t,x + f{t),y, ijj - 


- ft{t)y) 




{t, X, y, ip) ^ {t, x,y,'il) + g{t)), 





where G M and f{t) := e4,f{t) and g{t) := e^g{t). Accordingly, the admitted Lie symmetries of 
the barotropic vorticity equation on the beta-plane are scalings, time translations, translations 
in 2/-direction, generalized Galilean boosts in the x-direction and gaugings of the stream function 
with smooth time-dependent summands. 

We will compose transformations from these one-parameter Lie (pseudo)groups in the fol- 
lowing way r = Fg-^ o o T^g o Fj o T^ to a transformation F from the maximal Lie symmetry 
pseudogroup Gi of the vorticity equation (1). Any transformation of Gi then has the form 

(r,X,y,^) = (e^Hi + e2), e-'^{x + f{t)), e-^Hy + ^s), e-'^'^i^P + git)- ft{t)y)). (2) 

In what follows, we set h{t, y) = g{t) — ft{t)y for convenience and use the substitution hy = — ft, 
whenever hy occurs. 

Note that the maximal Lie invariance algebra go of the usual vorticity equation, which is also 
called the barotropic vorticity equation on the f-plane and corresponds to the value /3 = 0, is 
much wider than the algebra gi and contains gi as a proper subalgebra [1, 2]. The algebra go is 
spanned by the vector fields from gi jointly with the vector fields 

tdt - ipd^, -ydx + xdy, -tyd^ + txdy + + y^)c?^, h{t)dy + ht{t)xd^, 

where the parameter h runs through the set of smooth functions of t. This means that in addition 
to the transformations from Gi the maximal Lie symmetry pseudogroup Go of the usual vorticity 
equation also contains one more family of scalings, usual rotations in the {x, y)-plane, rotations 
depending on t with constant angle velocities and generalized Galilean boosts in the y-direction. 

The above form (1) of the vorticity equation is not particularly useful for a numerical evalu- 
ation. The reason is, of course, that any numerical model can be run only at a finite resolution. 



5 



which requires a suitably chosen averaging or filtering of Eq. (1). As from the point of view 
of invariant parameterization schemes the precise type of averaging is only of secondary impor- 
tance, we will employ a classical Reynolds averaging to Eq. (1) in the paper. This leads to the 
averaged vorticity equation 



where the right-hand side of this equation denotes the eddy-vorticity flux, which we aim to 
parameterize subsequently. For the sake of notational simplicity, we will omit bars over the 
mean quantities from now on. 

Slightly more generally, the vorticity equation (1) can be augmented with external forcings F 
and dissipative terms D yielding a general expression of the form 



A further question we aim to address is whether symmetries might be helpful in deriving invariant 
expressions for F and D. As by definition F denotes external forcing terms, it is not immediately 
clear why symmetries of the vorticity equation should place restrictions on the form of F. 
However, as we shall show, symmetries are valuable in finding invariant diffusion terms D that 
can be used in the course of turbulence modeling. For the sake of simplicity we therefore 
will use Eq. (4) for the case of F = and D 0, i.e. we assume that no external forcing 
acts on the system to which a damping is attached. Physically, the presence of F and D can 
be interpreted as symmetry breaking in the vorticity equation (1). Which symmetries are to 
be broken and which are to be preserved can be controlled upon expressing the term D via 
differential invariants derived in Section 5. This is a comprehensive problem and not all of the 
cases arising might be interesting from the physical point of view. We therefore restrict ourselves 
on the case where D or the eddy vorticity flux in Eq. (3) can be represented in such a manner 
that the resulting equation admits all the transformations from the maximal Lie invariance 
pseudogroup (2). This is the approach proposed in [32] and it appears to be suitable for the 
beta-plane equation. 

4 Algorithm for the construction of differential invariants 

Given a transformation pseudogroup G in the space of p independent variables x = (x^, . . . , x^) 
and q dependent variables u = (n^, . . . , u'^), in order to exhaustively describe its differential 
invariants one should find either a functional basis of differential invariants of any fixed or- 
der or a complete set of independent operators of invariant differentiation and a minimal set 
of differential invariants generating all differential invariants via invariant differentiation and 
functional combination. Within the framework of the method of moving frames the solution 
of this problem is split into two parts [14]. It is convenient to compute normalized differential 
invariants and operators of invariant differentiation using the explicit expressions for transfor- 
mations from G. The corresponding computation consists of two procedures, normalization and 
invariantization. At the same time, the derivation of syzygies (i.e., relations involving operators 
of invariant differentiation) between normalized differential invariants is mostly based on the 
determining equations of G, and an important tool for this is given by recurrence formulas. In 
this section we briefly describe related basic notions and results, paying the main attention to 
the computational realization of algorithms in fixed local coordinates. See [14, 19, 20, 36, 37, 38] 
for detail and rigorous presentations. 

In what follows the index j runs from 1 to p, the index a runs from 1 to q. We use two kinds of 
integer tuples for the indexing of objects. One of these kinds is given by the usual multi-indices 
of the form a = (ai, . . . , Up), where G No = NU {0} and |a| = ai + • • • + ap. By 6j we denote 
the p-index whose jth entry equals 1 and whose other entries are zero. Thus, both the derivative 




(3) 



Ct + tpxCy - ^yCx + P^x = F + D. 



(4) 
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9I"Im"/(5x^)"i • • • {dxP)"^ and the associated variable of the jet space J°°(x|n) are denoted by 
u%, D" = D"^ • • • Dp*", etc. Here Dj = d^j + Ylia a ^q+(5j'^«S operator of total differentiation 

with respect to the variable xK The other kind of index tuples is presented by J = {ji, ■ ■ ■ , Jk), 
where 1 ^ jk ^ p, k = 1,...,k, k € Nq. Such index tuples are used for the indexing of 
compositions of operators of invariant differentiation, which do not commute. Namely, we write 
Dj for D^-^ • • • Dj-^. The symbol dh denotes the horizontal differential, dh-F = Yl^=i(^j^)^^'' ^r 
a differential function F = F[u], i.e. a function of and u^. 

The normalization procedure for the pseudogroup G consists of three steps: 

1. Choose a parameterization (local coordinates) of G and find explicit formulas for the 
prolonged action of G in terms of the jet variables. 

2. Choose a subset of the transformed jet variables and equate the expressions for them to 
chosen constants. 

3. Solve the obtained system of normalization equations as a system of algebraic equations 
with respect to the parameters of the pseudogroup G including the derivatives of the 
functional parameters. 

The second step is nothing but a choice of an appropriate (coordinate) cross-section of the G- 
orbits. This should be implemented in a way ensuring that the system from the third step will 
be well defined and solvable. 

The normalization procedure results in the construction of a moving frame p for the pseu- 
dogroup G, which is, roughly speaking, an equivariant map from the jet space to G. Once the 
moving frame is constructed it can be used to map any object x(3^)^(n)) defined on an open 
subset of the jet space (a differential function, a differential operator or a differential form) to 
its invariant counterpart, t(x(ic, ""(n))) = xip^"'\x , u^^^)) " ix,U(^n)))- To carry out this in practice, 
one should replace all occurrences of the pseudogroup parameters in the transformed version of 
the object by their expressions obtained with the normalization procedure. 

Thus, the invariantization of the coordinate functions x^ and u'^ of the jet space yields the 
so-called normalized differential invariants = l{x^) and = i{u'^). In fact, the invari- 
antized coordinate functions whose transformed counterparts were used to set up the normal- 
ization equations are equal to the respective constants chosen in the course of normalization 
and hence these objects are called phantom differential invariants. Non-phantom normalized 
differential invariants are functionally independent and any differential invariant can be repre- 
sented as a function of normalized differential invariants. Invariantization of the operators of 
total differentiation, Dj, gives the operators of invariant differentiation, D^-, which upon act- 
ing on differential invariants produce other differential invariants. Note that the domain of 
the jet space, where invariantized objects are well defined, depends on what cross-section is 
chosen. 

In order to determine the algebra of differential invariants the normalized differential invari- 
ants and the operators of invariant differentiation play a key role. It has been proved [39] that 
for any Lie (pseudo)group the algebra of differential invariants can be completely described upon 
finding a finite generating set of differential invariants. As stated above, all the other differential 
invariants are then a suitable combination of the basis differential invariants or their invariant 
derivatives. The hardest part in describing the algebra of differential invariants is usually to find 
a minimal generating set of these invariants. Proving the minimality of a given basis usually 
involves the computation of the syzygies among the differential invariants, meaning functional 
relations among the differentiated differential invariants Dj-^^' ^i- ■ ■ '-^j-^o' . . . ) = 0. 

In general, the normalized differential invariants are derived from invariantization of the 
derivatives of the dependent variables, whereas the differentiated differential invariants are ob- 
tained by acting on normalized differential invariants of lower order with the operators of in- 
variant differentiation. The central point is that the operations of invariant differentiation and 
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invariantization of a differential function in general do not commute. Roughly speaking, the 
failure of commutation of these two operations is quantified by the so-called recurrence relations 

where co^ = L{dx^) [14, 37]. The forms ^-^ = and = t{(Pa) are the invariantizations of 
the coefficients of the general prolonged infinitesimal generator 

j=l a^Oa=l \ j=l / j=l 

of G. More rigorously, here ^-^ and are interpreted as coordinate functions on the space of 
prolonged infinitesimal generators of G, i.e., first-order differential forms in the jet space. Hence 
their invariantizations should also be forms, which are called invariantized Maurer-Cartan forms. 

The left-hand sides of the relations (5) are zero for phantom differential invariants. If the 
cross-section is chosen in a proper way, the recurrence relations for the phantom invariants can 
be solved for the independent invariantized Maurer-Cartan forms, which in turn can be plugged 
into the relations for the non-phantom differential invariants. Collecting coefficients of uj^ then 
yields a closed description of the relation between normalized and differentiated differential 
invariants, which in turn might enable the determination of a basis of differential invariants. For 
this latter task, specialized methods from computational algebra can be applied [38], which is, 
however, not necessary in the present case due to the relatively simple structure of the maximal 
Lie invariance pseudogroup Gi of Eq. (1). 

5 Differential invariants for the beta-plane vorticity equation 

In order to derive the moving frame for the maximal Lie invariance pseudogroup Gi of the 
barotropic vorticity equation on the beta-plane, it is necessary to prolong the group actions to 
the derivatives of ip. For this aim, we have to derive expressions for the implicit differentiation 
operators, Dy, Dx and Dy. They can be determined as the dual of the lifted horizontal coframe 
for Gi, which reads 

dhT = (Tf + iJtT^)dt + (T^. + ij^T^)dx + {Ty + ^PyT^)dy = e'^dt 

dhX = {Xt + iJtX^)dt + (X^ + ij^X^)dx + {Xy + ijyX^)dy = e~'^ ftdt - e"^^dx 

d^Y = {Yt + ^tY^)dt + (y, + ij^Y^)dx + {Yy + = e-"My. 

Therefore, the required implicit differentiation operators are 

Dr = e-^HD* - /tDx), Dx = e^^D,., Dy = e^^D^,, (6) 

where D^, D^,- and denote the usual operators of total differentiation with respect to t, x and y, 
respectively, = dt + Y,a^o,+s^d^^, = + Ea V'a+52 ^V'c and V)y = dy + Y.^'^c+Sid^^- Here 
and in what follows a = (01,02,03) is a multi-index running through Ng, |a| = ai + 02 + 03, 
5i = (1,0,0), 82 = (0,1,0), 82, = (0,0,1) and the variable il^a = 4'aia2a3 of the jet space 
corresponds to the derivative d^'^^ip / dt"^ dx°'^ dy°'^ . We also use the notation /(^^ = d^ f /dt^ and 
/i(^) = d^h/dt^, k G No. The transformed derivatives = 8^°'^'^/ /dT°'^dX°'^dY'^\ \a\ > 0, are 
then 

(D, - /,D.)-^o..«3 + j r^^"^^^^' = °' ' ]) ■ 

I ^(ai), a2 = as = j J 
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We carry out the normahzation procedure in the domain of the jet space which is defined by 
the condition ?/^^ 7^ 0. We choose the conditions for normahzation, which aUow us to express ah 
the pseudogroup parameters (including the derivatives of functional pseudogroup parameters) 
in terms of variables of the jet space: 

T = X = Y = 0, ^fcoo = ^feoi = 0, k = 0,l,..., ^010 = e, (7) 

where e = sgntpx- These conditions yield a complete moving frame for the maximal Lie invari- 
ance pseudogroup of the vorticity equation, which explicitly reads 

ei=lnv1V^, £2 = -t, £3 = -y, f = -X, 

= (Dt-VyD^)^V'3/, h^k) = -{^t-^y^xr^, /C = 0,1,.... 

The series of equalities for f[k+i) ^-^i ^(fc) is proved by induction with respect to k using the 
equations 

/(fc+i) = (Dt - ft-D^f^y, h^k) = -(Dt - /tD,)V. 

The nontrivial normalized differential invariants are found via invariantizing the deriva- 
tives ipa for the values of a for which are not involved in the construction of the above 
moving frame, i.e., for 

a G A = N|] \ {{k, 0, 0), [k, 0, 1), (0, 1, 0), G No}. 

In other words, for each a € ^4 we should substitute the expressions (8) for the pseudogroup 
parameters into the expressions for (The invariantization of the coordinate functions chosen 
for the normalization conditions (7) are equal to the corresponding normalization constants and 
are the phantom normalized differential invariants for the moving frame (8).) As a result, we 
obtain the differential invariants 

la = iii'a) = |Vx|("^+"^-"^-')/'(Df - ^l^yD^T'^Oa.a,. Oi ^ A. 

The order of Iq as a differential function of "0 equals |q;|. It is also obvious that any finite number 
of the invariants la are functionally independent. This agrees with the general theory of moving 
frames [14, 19, 37], which also implies a stronger assertion. 

Theorem 1. For each r ^ 2 the functions 1^ = \ipx\^"''^~^°'^~°'^^^^^'^(Dt — ipy^x)°'^'4^0a2a3, where 
a £ A and \a\ ^ r, form a local functional basis of differential invariants of order not greater 
than r for the maximal Lie invariance pseudogroup Gi of the barotropic vorticity equation on 
the beta-plane. 

The description of differential invariants of Gi given in Theorem 1 is sufficient for applica- 
tions within the framework of invariant parameterization. At the same time, it is interesting and 
useful to have more information on the structure of the algebra of differential invariants of the 
pseudogroup Gi including the operators of invariant differentiation. A complete set of such in- 
dependent operators is derived by invariantization of the usual operators of total differentiation, 
yielding 

Dj = ^l=(Di-V,D,), DL = ^^D,,, Dl = y^\Dy. (9) 

This is practically realized via substituting the expressions (8) for the pseudogroup parameters 
into the implicit differentiation operators (6). Any operator of invariant differentiation related to 
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the pseudogroup Gi is locally a combination of the operators (9) with functional coefficients de- 
pending only on differential invariants of Gi. The commutation relations between the operators 
D\, and Dy are 

[Dj,Dy = |lo2oDi + (/on + |/iio) D^, 

[Dj, = |/oiiDj + /002DL + |/iioD;^, (10) 

[DL,Dy] = -I020DJ; - 2^oiiDL- 

In order to completely describe the algebra of differential invariants of Gi , it remains to estab- 
lish a basis of differential invariants such that any differential invariant of Gi can be represented 
as a function of basis elements and their invariant derivatives. It is also necessary to compute a 
complete system of syzygies between basis invariants. For this aim, we will evaluate the recur- 
rence relations between the normalized differential invariants and the differentiated differential 
invariants as detailed in [14, 37]. The starting point for the application of the general algorithm 
to the maximal Lie invariance pseudogroup Gi of the vorticity equation on the beta-plane is the 
system of determining equations for the coefficients of a vector field 

Q = T{t, X, y, i))dt + i{t, X, y, ip)d:c + r]{t, x, y, il;)dy + ip{t, x, y, ip)d^ 

from the maximal Lie invariance algebra of Eq. (1), which reads 

Tx = Ty = = = = rjt = = Tj^ = (fx = 0, 

ix='ny = -Tt, (Py = -£,t, (Pr(, = -3rf. 

Consider the prolonged operator Q^o = Q + Yl\a\>o ^"'^i^a- The coefficients of Qoo are calculated 
by the standard prolongation formula [35, 39], 

In view of the determining equations, the coefficients ip" take the form 

^-=(a.+.,-a,-3)n*<.-|:(';')f'"*"-"--'"+{ ;,tr" Z=l=o'}- 

where ^(fc) = d'^S^/dt'^ and (/^(fc) = d^(p/dt^, k = 0,1,2,.... We collect the coefficients of Q 
and their derivatives appearing in the expressions for the prolonged coefficients of Q and de- 
note the associated invariantized objects, which are differential forms, as f'' = i-(r), = /.(rt), 
= i{^[k))i V = '•(^) and (^'^ = i[ip(^f.-^). In the course of the normalization (7) the invariantized 
counterparts = i{(p°') of the prolonged coefficients of Q are 

k=i ^ ^ k=i ^ ^ 

0" = ("2 + as - ai - 3)Iaf^ ~ X] ( k) ^a-kSi+Sii'' if (22 > or > 1. 

k=i ^ ^ 

For lower values of |a|, < |a| ^ 3, we calculate 



^100 




(^010 = 


-2f\ 


^ = -4 , 




^200 


=0'- ee 






= -3/lloT"'^ - Io20^^, 


-101 t2 T n 

= -? - -'OIK 


^020 




^011 = 




-002 T ~1 
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^300 


= ^'- 


- 3/iiof - 




= -4/2ior^ 


- -^020^^ - 




= -2/i2ot' 






-021 -012 -003 


^030 


= =ip 


= if -- 



-201 



-111 or ^1 r £1 ,?,102 or r tl 

95 = -2iiiir - io21? , ^ = -^ho2T - i012? , 

= 0. 

From the recurrence relations for the phantom invariants = = i{x), = L{y), 

Im = L-{ipm), Im = L-{ipm), i = 0, 1, . . . , and /qio = i(V'oio), 

dh//° = + f° = 0, dh//^ = + e° = 0, di,H^ = oj^ + fi = 0, dh/ooo = + = 0, 



dh/joo = IjioUJ^ + ^ - ee - (^j] Ij-k,wi^ = 0, j = 1, 2, . . . , 

k=i ^ 

dh/ioi = /iii^' + /i02a;=^ - - ^ /j-fc,iif = 0, j = 0, 1, . . . 



k=l 

dh-^010 = -Z'lio^"^ + -^020^2 + loiiuj^ - 2f^ = 0, 
we derive expressions for the invariantized Maurer-Cartan forms 



r 



e = + /j-i,02w^ -^(\ 

k=l ^ ^ 



k=l 

j = 1,2, . . . . The forms should be calculated recursively starting from j = 1. Thus, 

= -^011"^^ + -^002"^^, 

? = (All — -^oii)^^ + {ho2 - -^011-^002)^^, 

= (-^211 — 3/011/111 + /fii)i^2 _|_ (/2Q2 — 3/011/102 + /oil-^002)'^^, 

In general, ^-^ = ^^^'^uP' + ^■^'^w'^, where the coefficients ^■''^ and ^■''"^ are expressed in terms of 
normalized invariants /^ with \a\ ^ j ' + 1. 

The recurrence relations for non-phantom normalized invariants correspondingly read 

- ^ f ^ j la-kSi+Sii'' if 02 > or as > 1. 
k=l ^ ^ 

As by definition d^F = {D\F)uj^ + (D!j_.F)a;2 + (DJ^/*")!:!;^, the above recurrence relations can 
be split into a list of equations for first-order invariant derivatives of normalized differential 
invariants /„ with 02 > or 03 > 1 by taking into account the expressions for the invariantized 
Maurer-Cartan forms: 

r>i r - 7- i "2 + Q3 - ai - 3 

>->t^a — ^a+5i H ^ -'llO-'aj 

ni r - r j. «2 + Q3 - Qi - 3 V" /^"i V /?fc.2 

■L'x-'a — -'q+52 "I 2 ^020^a - ^ JJ^a-k5i+S2Z > ^-^-^^ 

r>i r - r , «2 + Q3 - Qi - 3 ^ /^"A r ;?fc.3 

■Uj^-'a — -'a+53 H 2 -'Oll-'a ~ 2^ I ^ l-'a-fc5i+52? • 

z — 1 V / 
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We only present the closed expressions for the first-order invariant derivatives of Iq with \a\ ^ 3: 



Who 


— -^210 


3 r2 p 
2-'ll0' ^ 


'L-^llO — /12O — |/llo/o20 — /0I1/02O; 






Dy/llO 


= /lll 


— f -^110-^011 


- /020/002, 






Dt/o20 


= -^120 


— ^-^110-^020, 


D^/020 — /030 — ^-^0201 Dj^/020 — /02I — |/01l/020) 


Dj/oii 


= /lll 


— ^-^110-^011 1 


D^/oii = /021 — ^/oii/o20i Dy/oii = 


- /012 — 


2-'011i 


Dj/002 


= -^102 


— ^-^110-^002, 


D^./002 = /012 - ^/020/002, D'y/002 = 


= /003 — 


^/01l/002, 


D(/210 


= -^310 


— 2/110/2IO) 


D^/210 = /22O — 2/020/210 — 2/011/120 


+ (4i 


— /lll)/020 


DJ,/210 


= -^211 


— 2/011/210 ■ 


- 2/002/120 + (/002/0II — /l02)/020; 






0^/201 


= -^301 


— 2/110/201 ) 


D^/201 = /2II — 2/020/201 — 3/001/111 


^ -'oil! 




DJ,/201 


= -^202 


— 2/011/201 ■ 


- 2/002/111 — /0I1/102 + /o02/oil5 






Dj/l20 


= -^220 


— /1I0/12O) 


D|j./l20 = /l30 — /020/12O — /01l/030 5 






Dj,/i20 


= -^121 


— /0I1/12O - 


/002/0305 








= -^211 


— /iio/iii> 


DL-^111 = /121 — /020/111 — /011/021, 






D'/iii 


= A 12 


— /oii/iii - 


/002/02I, 






0^/102 


= -^202 


— /110/102, 


D|j./io2 = /112 — /020/102 — /011/012, 






D^/102 


= ho3 


— /oi 1/102 - 


/002/012, 






Dj/oso 


= -^130, 


D^/030 = 


/040> D^/030 = /031; 






Dj/021 


= -^121, 


D^/021 = 


/031i Dy/021 = /022, 






Dj/012 


= -^112, 


D!j./oi2 = 


/022, D'y/012 = /013, 






Dj/o03 


= -^103, 


D!j./o03 = 


/013i D'y/003 = /004- 







In principle, it is possible to read off the generating differential invariants from the above split 
recurrence relations. The expressions for la+Sn /«+52 ^a+h derived from (11) only involve 
first-order invariant derivatives of /q and normalized invariants of orders not greater than |a|. 
This implies that a generating set of differential invariants consists of invariantized derivatives 
which are minimal with respect to the usual partial ordering of derivatives and are not phantom 
invariants. We have four such minimal elements, 

-/no — , — , Im.n — — , -irm — — , -fO"~ 



ViV'xP ' ^A^\ 

All the other invariantized derivatives are expressed via invariant derivatives of /no, /o205 /on 
and /oo2- As was indicated above, not all differentiated differential invariants are necessarily 
functionally independent, which is encoded in syzygies of the algebra of differential invariants. 
Taking into account these syzygies can further reduce the number of generating differential 
invariants thereby allowing one a more concise description of the basis of differential invariants. 
In the present case, we find the following lower-order syzygies: 

Df/oii — Djy/110 = /no/on + /020/002, 
Dj/020 — D!j,/iio = /o2o(/iio + /on), 
Dy/011 - Di;/oo2 = ^/o2o/oo2 — ^/qh, 
DJ^Joii — Dy/020 = 0, 

(Dy)^/iio — DjD^/002 = ^(Df — /oii)(/o2o/oo2) - (DJ/ + /oii)(|/iio/oii + /020/002) 

— /oilDy/110 — /002DJ//020, 
(Dy)^/020 - (D!c)^-^002 = ^D!e(/020/002) — ^Dj/(/oil/o2o)- 
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From the two first syzygies we can express the invariants Iqii and /002 via invariant derivatives 
of /no and /020, 



-'oil — J -'110) 

-'020 



1 .^.i , / DU020-DU110 , \ d;,Jiio 

l^f - -'llOj r -'110 r 



-^002 — 7 (Df — -^110) 1 r '^'^^ I T 

^020 \ ^020 J -'020 

Another way of finding relations between generating invariants is to use the commutation 
relations between the operators of invariant differentiation. Evaluating each equality from (10) 
on an element I from the above generating set, we obtain a system of linear algebraic equations 
with respect to the other elements of these sets, which can be solved on the domain of the jet 
space where the determinant of the matrix associated with the system does not vanish. It is 
convenient to choose, e.g., / = /o20- Then, we derive the representations 



-^02oDy-fo20 — 2e[D^, D^]Io2o 

Ux-'020 

2e[DJ,D^]/o2o - -fo2oDj/o20 
-/no — ;:rr-:^ ^£-(0 



'110 — T^fT ^t-'Oll, 

U^-/020 

_ [Di,Dyio2o eDj/020, eDj^/020, 

U^i020 ^ Ux-'020 ^ -L'^-'020 



which are defined on the domain of the jet space where D|j./o20 7^ 0, i.e., D^(-y/|i/IjJ| ) 7^ 0. 
As a result, it is straightforward to establish the following theorem. 

Theorem 2. The algebra of differential invariants of the maximal Lie invariance pseudogroup 
of the barotropic vorticity equation on the beta-plane (1) is generated, in the domain Qi of the 
jet space where Dj(y^|^/j^) ^ 0, by the single differential invariant /020 = V'xx-/ \/|V'x-| along with 
the three operators of invariant differentiation 

1 



All other differential invariants are functions of /020 and invariant derivatives thereof. 



6 Invariantization of parameterization schemes 

The Replacement Theorem states that any differential invariant /(x,it(„)) of order n can be ex- 
pressed in terms of the normalized differential invariants via replacing any argument of I{x, U(„)) 
by its respective invariantization, see [20]. In particular, any system of differential equations can 
be represented using the normalized differential invariants of its associated maximal Lie invari- 
ance group. The invariantization of the vorticity equation (1) in view of the moving frame (8) 
reads (I120 + ho2) + (-^021 + -^003) + /? = 0, or, explicitly 

^l^M^ + Cy + f3 = 0. (12) 

This is the fully invariant representation of the barotropic vorticity equation on the beta-plane. 

Differential invariants computed in the previous section can be assembled together to invariant 
parameterizations of the eddy- vorticity fiux in the averaged vorticity equation (3). Alternatively, 
we can invariantize any existing parameterization scheme under the moving frame action (8). 
The following two examples implement this idea. 
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Example 1. A classical albeit simple parameterization for the eddy-vorticity flux is 



evf := {C%)^ - (C'V;), = D,.(KC.) + D,(KC,), 

where K = K{x, y) might be considered as a spatially dependent function. The most straightfor- 
ward way to cast this parameterization into the related invariant one is by applying the moving 
frame (8) to the terms on the right-hand side. This yields 

evf = D!,.(i^(/o30 + /012)) + D|^(i^(/o2i + /003)) = K{Io4o + 2/022 + /004) 

= Ky^llp^KCxx + Cyy), 

where evf = i(evf) and K = const now as l{x) = L{y) = 0. The invariant representation of the 
closed barotropic vorticity equation then reads 

^i^Ml + Cy + f3 = Ky^liCxx + Cyy). 

Example 2. The anticipated (potential) vorticity method was originally proposed by Sadourny 
and Basdevant [46] . The idea of this method is to approximate the diffusion effect in the vorticity 
equation by a weighted upwind estimate of the vorticity itself, i.e. by employing 

7]t + J{iP,r]) =i/J(V^,A"J(V',r/)), 

where is a constant, n € Nq and 7] is the absolute vorticity. Here and in what follows A = 
is the two-dimensional Laplacian. The purpose of adding the specific forcing term on the right- 
hand side of the vorticity equation is to suppress the high frequency noise in the vorticity field 
and at the same time to ensure that energy is conserved during the integration while enstrophy 
is dissipated. The latter properties can be easily verified upon multiplying Eq. (1) with the 
stream function ip and any function of the absolute vorticity 77, respectively, and integrating 
over the domain fi, see also [54]. 

There is a problem with this parameterization scheme in that it is not Galilean invariant. 
Galilean invariance (as well as the proper scale invariance), however, can be easily included by 
the method of invariantization. For the sake of demonstration, we consider the case of n = 
here, which is the original version of the anticipated vorticity closure. Upon using the moving 
frame (8), we obtain 



L{J{ij,J{Tp,r]))) = —==J{'4)y,ri) + VW^\r]yy. 

Attaching this to the invariant representation of the vorticity equation (12), the vorticity equa- 
tion with fully invariant closure reads (e = sgnipx) 



Vt + Jii^,v) = vyWAieJii^y^r]) + i)xVyy)- (13) 

It is obvious that this parameterization is quite different from that proposed in [46]. It cannot 
be brought in the form of nested Jacobian operators and it does not conserve energy any more 
(for the derivation of conservative invariant closure schemes, see Section 8). On the other hand, 
the inherent invariance of the closed vorticity equation (13) under Galilean and scale symmetry 
is an appealing property for itself and might be relevant e.g. when vorticity dynamics is studied 
in a moving coordinate frame. 

Quite recently, an approximate scale invariant formulation of the anticipated potential vor- 
ticity method was proposed in [15] using scale analysis techniques and physical reasoning. The 
motivation for this study is that modern weather and climate models might be required to 
operate on grids with variable resolution. Unfortunately, varying resolution in an atmospheric 
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numerical model is not a simple task as most of the parameterization schemes employed are def- 
initely not scale invariant, but rather tuned to yield best results on some fixed grid. This means 
that painful efforts might be necessary in order to adjust all the parameterization schemes of 
a numerical model to different spatial-temporal resolutions. Having a general method at hand 
that allows deriving of scale insensitive closure schemes is therefore of potential practical interest 
in numerical geophysical fluid dynamics. Albeit simple, the method of invariantization of ex- 
isting parameterization schemes may give appropriate closure schemes that are both physically 
meaningful and respect those symmetries that might be essential for a specific process to be 
represented numerically. 

These are only two examples for fully invariant closure schemes. See one more example in the 
next section. In principle, each term of the form S{I^, . . . ,1^), where S" is a smooth function 
of its arguments and I^, . . . , are differential invariants of Gi, satisfies the same requirement 
when added to the right hand side of Eq. (12). In other words, the general form of closure 
ansatzes for Eq. (12), which are invariant with respect to the entire group G\, is 

7 Application of invariant parameterizations to 
turbulence modeling 

In this section, we give an application in which we aim to demonstrate in practice the ideas 
outlined above and in [41]. This example deals with turbulence properties of the two-dimensional 
incompressible Euler equations. Strictly speaking, turbulence is a three-dimensional problem as a 
two-dimensional turbulent flow is not stable with respect to fully three-dimensional perturbations 
to that flow [47]. Nevertheless, there are countless studies concerning the turbulent properties 
of two-dimensional flow simply because it is a relevant problem in large-scale geophysical fluid 
dynamics, which behaves as approximately two-dimensional. 

In short, the first theoretical results concerning two-dimensional turbulence were derived in [4, 
28], following the pioneering work on three-dimensional turbulence done by Kolmogorov [27]. 
Extensive numerical studies have been carried out since then attempting to verify distinct aspects 
of the theory proposed [6, 7, 13, 21, 30]. The two-dimensional case is especially peculiar, as it 
admits infinitely many conservation laws including the conservation of energy. The energy in 
the barotropic vorticity is purely kinetic and can be represented in different ways using doubly 
periodic boundary conditions as 



where il. = [0,Lx[x [0,Ly[ and d^ = dxdy. The special form of Eq. (1) leads to the following 
class of conservation laws 



for any smooth function g of the absolute vorticity r] = ^ + fo + /32/. The most relevant realization 
of the above conservation laws in the present context is the enstrophy, given for the particular 
value g = if' /2. 

First of all, consider the case of no differential rotation (/3 = 0), i.e. the Coriolis parameter f 
is approximated by the constant fo, which is referred to as the f-plane approximation. It is the 
simultaneous conservation of energy and enstrophy in this case that leads to the remarkable 
behavior of two-dimensional turbulence [47, 53]. Starting with a random initial velocity (or 
stream function field), energy is transported to the large scale, while enstrophy is transported to 




(14) 
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the smaller scales. This cascade is associated with an organization of the vortices, with vortices of 
the same sign merging into bigger ones (though the precise mechanisms of the cascade including 
the role of the vortices are not yet fully understood). In order to initiate these fluxes of energy 
to the larger scale and enstrophy to the smaller scale and thus the process of organization, it 
is necessary to place a sink of enstrophy at the very small scales. This sink acts as a remover 
of enstrophy while ideally conserving energy, as the latter is transported away from the small 
scales on which the dissipation acts (which in practice is hard to realize in a numerical simulation 
using a finite number of grid points). It is believed that the form of the energy spectrum in a 
range above which dissipation is acting (inertial range) can be derived using scaling theory in a 
similar manner as it was shown by Kolmogorov for the three-dimensional case. 
The energy spectrum E{k) is defined by 

£ = -i- / v^d^ = — i- / {V^l^fdA = [ E{k)dk, 

where £ is the average energy, k = \J (A;^)^ + (/cJ/)^ is the scalar wave number, k^ and k^ are the 
wave numbers in x- and y-direction, respectively. The possibility of using a single wave number 
is due to the assumption of isotropy that is generally made in turbulence theory and which is 
reasonable in the case of vanishing differential rotation. According to the theory, the form of 
the energy spectrum in the inertial range should follow 

E(k) oc 

This is referred to as the enstrophy cascade in two-dimensional turbulence. 

The impact of the beta-term in the vorticity equation on the turbulent cascades was first 
studied in [45]. In this seminal paper, it was remarked that the Rossby wave solutions admitted 
by the beta-plane equation can act as a source of anisotropization of turbulence at the larger 
scale. Qualitatively, at some stage the size of the vortices is big enough that they are exposed 
to the effect of differential rotation, which essentially hinders the tendency of vortex growth due 
to the inverse energy cascade. Rather, the vortices evolve into Rossby wave and eventually to 
the formation of zonal jets as observed e.g. on giant planets. Depending on the precise setting 
used (e.g. strength of the differential rotation, additional energy injection to the system), the 
results of turbulence simulations can vary, but often energy spectra steeper than those predicted 
theoretically can be found [23, 31, 45]. 

In practice, the sink of enstrophy at the small scales is usually implemented by adding a 
hyperviscosity of the form 

D = (-1)"-VA"C (15) 

for n G N"*" to the right-hand side of Eq. (1), cf. Eq. (4). However, it can easily be checked 
that this form of hyperviscosity is not invariant under the Lie symmetry pseudogroups of the 
beta-plane and f-plane equations. More specifically, it violates the scale invariance of Eq. (1). 
From the theoretical point of view, this violation appears to be especially odd, as it is precisely 
the scale invariance of the Euler equations that is used to derive the form of the energy spectrum 
in the inertial range. 

Theorem 1 directly implies that the invariantization l{D) = (— l)"~^i/y^]^0^^5p*^A"C is a 
differential invariant of the maximal Lie invariance pseudogroup of the vorticity equation. In 
view of the results of Section 6, we conclude that the form of the diffusion term obtained in the 
course of the invariantization is 

b = = (-i)"-vviV'x|2"+iA\. 

The completely invariant formulation of the vorticity equation on the beta-plane with hyper- 
diffusion therefore reads 

Ct + ^xCy - ^yCx + P^x = (-l)"-VVlV'x|2-+lA\. (16) 
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Note, however, that the price for introducing an invariant enstrophy sink is the nonlinearity 
of the (hyper)diffusion term. More generally, the situation is alike to the problem of finding a 
relation between the Reynolds stresses and the mean strain rate in the Reynolds averaged Navier- 
Stokes equations or in large-eddy simulations thereof. It was pointed out that establishing a 
relationship between the nonlinear Reynolds stresses and the linear strain rate (i.e. invoking 
the Boussinesq hypothesis) may lead to unrealistic results for certain turbulent flows such as 
in rotating or stratified fluids or those exposed to abrupt changes of the mean strain rate, see 
the discussions in [40, 55]. It is therefore worthwhile pointing out that the requirement of 
preserving the entire maximal Lie invariance pseudogroup of the barotropic vorticity equation 
on the beta-plane automatically yields nonlinear hyper diffusion terms. For n = 1, the right- 
hand side of Eq. (16) can be considered as a generalized down-gradient parameterization for 
the eddy- vorticity flux, which is also a nonlinear quantity. That is, requiring a (hyper) diffusion 
scheme to be scale invariant, it is indispensable to use nonlinear (hyper)diffusion. 

It is important to note that the anisotropic coefficient Y^|?/^^p"+^ arises due to the special 
form of normalization conditions (7) we have chosen in Section 5 for the construction of the 
moving frame. This form is by no means unique but rather a consequence of the moving frame 
we have invoked. The situation is comparable to the discretization of differential equations, 
which can also be done in multiple ways. Some schemes have better properties than others 
and ultimately it is necessary to both analyze and test the various schemes for different sets of 
problems. Having more than one possibility to construct invariant subgrid-scale schemes out 
of a given non-invariant scheme should therefore be considered as an advantage rather than 
as a drawback of the proposed method. The knowledge of the complete set of differential 
invariants, which is obtained as a byproduct when determining the invariantization map for a 
given group action, allows one to derive series of invariant closure schemes starting from that 
obtained as a direct result of the invariantization of the given initial scheme. This is facilitated 
by recombining a given invariant scheme using the differential invariants, as any functional 
combination of differential invariants is again a differential invariant. 

A number of alternative (isotropic) forms of a completely invariant nonlinear hyperviscosity 
term for the vorticity equation on the beta-plane can therefore be suggested, e.g. 

i) = (-l)"-ii.C2n,+i^n^^ 2) = (-l)"-VV(C^"+iVA"~iC), etc., 

which are derived upon recombining the differential invariants derived in Theorem 1. Due to 
the wide possibility for varying ansatzes for invariant parameterizations we can control differ- 
ent desirable conditions which proper invariant closure schemes should additionally satisfy, cf. 
Section 8. 

Subsequently we will exclusively work with Eq. (16). Our motivation for choosing the 
anisotropic hyperdiffusion (16) rather than any of the above isotropic ones stems from recent 
experiments on turbulence which suggest that contrary to the Kolmogorov hypothesis the small 
scales might indeed feel the effects from the large scale being anisotropic, i.e. that anisotropy can 
propagate through to the very small scales, see e.g. [49]. However, future tests will be conducted 
so as to compare the different forms of invariant hyperdiffusion. 

Remark 1. In order to set up a numerical model, a decision has to be taken about which 
boundary conditions should be implemented. It is very rare that the numbers of symmetries 
admitted by a differential equation is not reduced for an associated boundary value problem. The 
most immediate boundary conditions in the atmospheric sciences are periodic ones. However, a 
periodic domain implies a fixed domain size and therefore breaks the scale invariance of Eq. (1). 
On the other hand, scale invariance is an equivalence transformation of the class of all periodic 
boundary value problems of the vorticity equation, see also the discussion in [11]. A more serious 
problem is that the periodicity in y-direction is not natural for the beta-plane. On the other 
hand using a channel model (rigid walls in the North and in the South of the domain) breaks also 
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Figure 1: Energy spectra derived upon integrating Eq. (16) for n = 2, normalized by the 
respective values at wavenumber fc = 1 for the sake of comparison. 



the translational invariance in y-direction thereby reducing the admitted Lie symmetry group 
even stronger than in the presence of doubly periodic boundary conditions (though, in contrast 
to usual hyper diffusion, it would not be necessary to define an additional boundary condition for 
the invariant hyperdiffusion as by definition -i/^^; = at the walls of the channel and the diffusion 
term therefore vanishes there). This is why we will use doubly periodic boundary conditions 
although /3 7^ here. Despite this slight inconsistency, doubly periodic boundary conditions are 
used quite extensively in studying turbulence properties on the beta-plane [31, 45, 53]. 

We give some numerical experiments using Eq. (16) and integrating it with a finite difference 
scheme. An invariant biharmonic dissipation is used in all the experiments, i.e. n = 2. The 
nonlinear terms on the left-hand side are discretized using the Arakawa Jacobian operator [3], 
which guarantees energy and enstrophy conservation of the spatial discretization in the case of 
vanishing dissipation, z/ = 0. A leapfrog scheme is used for the time stepping in conjunction 
with a Robert-Asselin-Williams filter [56], in order to suppress the computational mode. The 
size of the domain is = Ly = 2.56 • 10^, with a default of = 128 grid points in each 
direction, j3 = 1.6 • 10~^. The initial condition is a Gaussian random stream function field, with 
the initial energy spectrum given by the function E{k) oc /{I + /e/32)^^. The physical setting 
is equivalent to shrinking the Earth to approximately one percent of its original size. 

In Fig. 1 the scale invariance of the energy spectra is shown after 4320 time steps of integra- 
tion. Scale invariance is verified by scaling t, x, y and ip using 

Pe, : {t, X, y, i;) ^ {e^H, e'^'x, e-'^y, e'^'^^p), 

for the the values ei = 1, ei = 2 and £i = 1/2. This invariance is of course to be expected from 
the invariant formulation of Eq. (16). The difference between the single runs is due to numerical 
errors. The resulting stream function fields also coincide (not shown here). 

The corresponding integration with the noninvariant hyperdiffusion scheme is shown in Fig. 2 
using scale factors ei = 1 and ei = 2 (the corresponding integration for ei = 1/2 already 
becomes unstable in the setting used). Clearly, the spectra are not invariant (nor coincide the 
stream function fields), which is due to breaking of the geometric structure in the vorticity 
equation. We argue that this is not natural as the scale invariance of the incompressible Euler 
equations should imply the lack of any fixed scale, which is not realizable upon using usual 
hyperdiffusion. 
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Figure 2: Energy spectra of Eq. (1) upon using the usual hyperdiffusion (15) for n = 2, normal- 
ized by the respective values at wavenumber k = 1 for the sake of comparison. 

Besides invariance, another relevant question in the present framework is whether the form of 
the spectra derived is reasonable from a physical point of view. This question is of course much 
harder to answer as there are several open problems in two-dimensional turbulence such as the 
precise nature of the interaction of the vortices among each other or with the background field. 
Nevertheless, we shall like to present the energy spectrum for fully developed freely decaying 
turbulence using Eq. (16). To this end, we conduct a numerical experiment with N = 512 grid 
points in each direction and integrate the equation over approximately 10000 time steps. The 
result of this integration is shown in Fig. 3. Qualitatively, attaching the invariant hyperdiffusion 
to the vorticity equation leads to an increase of the energy in the begin of the integration (in our 
simulation the energy initially grew about 16%). After the initial growth the system stabilizes 
and after a while starts to dissipate energy (and enstrophy). The spectrum shown in Fig. 3 was 
sampled approximately at the onset of the stabilization of the energy level after the transient 
growth phase. Note, however, that the nonlinearity of the diffusion term makes it somewhat 
hard to evince the universality of the derived spectrum, which becomes even harder than in the 
case of usual decaying turbulence [23]. We should like to stress, though, that the spectra we 
obtained by performing different runs all resemble the one presented in Fig. 3. The lower part 
of the spectrum (above the dissipative range) follows fairly good the theoretically predicted 
form, while the upper part is slightly shallower with slope of approximately —2.7. Shifting the 
initial energy spectrum E{k) towards smaller wave numbers leads to spectra that display a —3 
slope over a wider range of wave numbers in our simulations (not shown here). This points to 
the well-known behavior that the form of the inertial range may depend on the initial conditions, 
see also the discussion in [5, 48]. 

Fig. 4 shows the associated stream function at the end of the integration. When vortices 
grow in size, they start to feel the effect of differential rotation and are elongated in East- West 
direction. As is well known, the elongation causes anisotropization of turbulence at the larger 
scale. The anisotropization is a further reason why spectra of E{k) against k should be taken 
with care after sufficiently long time of integration as the assumption of isotropy justifying the 
usage of a single wave number might no longer be valid. 

Remark 2. Decaying turbulence simulations are an important class of tests for numerical inte- 
gration schemes. On the other hand, from the point of view of both the theory and application, 
it is generally more instructive when Eq. (1) is augmented with some forcing which supplies 
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Figure 3: Energy spectrum derived upon integrating Eq. (16) for n = 2. The solid line denotes 
the k~'^ spectrum predicted by the theory for the enstrophy inertial range, the dashed line 
follows A;"^-'^. 




Figure 4: Final stream function after approximately 10000 time steps of numerical integration 
of Eq. (16) for the decaying turbulence experiment. The elongation of large-scale vortices in 
East-West direction due to rotation is clearly visible. 
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energy to the system and thereby prevents turbulence from dying out. As it is then usually 
necessary to damp out the energy which is otherwise piling up at small wave numbers (large 
scales) due to the inverse energy cascade, an additional drag term is introduced in Eq. (1). This 
drag term can be either physical (e.g. linear Ekman drag due to bottom friction) or, similar as 
hyperviscosity, scale selective. In the latter case, one uses a hypoviscosity [16], which is given 
by adding a term proportional to A~"^, which acts scale selective by emphasizing the large 
scale and thus is effectively energy removing. Again, one could raise the question whether such 
a hypofriction should possess some invariance properties, but this is beyond the scope of the 
present paper and should be considered in a forthcoming study. 

8 Conservative invariant parameterizations 

A parameterization is called conservative if the corresponding closed system of differential equa- 
tions possesses nonzero conservation laws. Special attention should be paid to parameterizations 
possessing conservation laws that have a clear physical interpretation (such as the conservation 
of energy, mass, momentum, etc.) and that originate from the conservation laws of the initial 
system of equations. If a parameterization is both conservative and invariant with respect to a 
pseudogroup of transformations, it is called a conservative invariant parameterization. 

The general method for singling out conservative parameterizations among invariant closure 
ansatzes is based on the usage of the Euler operators, i.e. variational derivatives with respect to 
the dependent variables [35]. Suppose that Cg: U{x, 0) = 0, I = I, . . . ,m, 9 = 6{I^, . . . , I^) 
represent a family of local parameterizations for a system C: L\x, U(„)) = 0, / = 1, . . . , m, which 
are invariant with respect to a pseudogroup G. Here V are fixed smooth functions of their 
arguments. The tuple 9 of arbitrary elements consists of smooth functions of certain differential 
invariants I^, . . . , of G. It runs through a set of such tuples constrained by a system of 
differential equations, where I^, . . . , play the role of independent variables. We require the 
tuples (A*"^, . . . , A"''), m = 1, . . . , M, of differential functions of u to be characteristics of M 
linearly independent local conservation laws of the system Cq for some values of 9, i.e. for each m 
the combination X^^L^ + • • • + A™'L' is a total divergence. The theorem on characterization of 
total divergences [35, Theorem 4.7] then implies that 



for each m = 1, . . . ,M and a = 1, . . . , g, where is the Euler operator associated with the 
dependent variable n", E"/ = X]cj(— D)"/u;^- Splitting Eqs. (17) with respect to derivatives of u 
wherever this is possible, one constructs the system of determining equations with respect to 9, 
which should be solved in order to derive the corresponding conservative invariant parameteri- 
zations. 

As the direct computation is too cumbersome, we use some heuristic arguments and look 
for a diffusion ansatz for the barotropic vorticity equation on the beta-plane which satisfies the 
following relevant and valuable conditions: 

• It is invariant with respect to the entire maximal Lie invariance pseudogroup Gi of Eq. (1). 

• The subgrid-scale term or, more generally, the sink term to be represented is a differential 
function of the vorticity. 

• This expression is as similar as possible to the hyperviscosity term (15). 

• And, last but not least, the parameterization is conservative. More precisely, it possesses 
all the conservation laws of Eq. (1) with zero-order characteristics. 

An example of such a parameterization is given by 



(17) 




(18) 
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All the properties listed above can be checked for the sink term (18). Thus, the expression 
for D from (18) involves only the vorticity and its derivatives and is quite similar to (15). 
Moreover, the diffusion is a globally defined differential function which is a polynomial of 
its arguments. The invariance of Eq. (18) with respect to Gi can be simply checked using 
the infinitesimal invariance criterion. A more sophisticated way to check this invariance is to 
rewrite Eq. (18) in terms of normalized invariants of the pseudogroup Gi, which will not be done 
explicitly here. As an unexpected but valuable bonus we have that the maximal Lie symmetry 
pseudogroup of Eq. (18) with the same term D in the case of the f-plane {j3 = 0) is even wider 
than Gi. It also includes the usual rotations of the variables (x,y) and the generalized Galilean 
boosts in y-direction, which belong to the Lie symmetry pseudogroup Go of the barotropic 
vorticity equation on the f-plane. This in particular means that the parameterization (18) is 
isotropic. 

The space of zero-order characteristics of Eq. (1) is generated by the characteristics A = /(t), 
A = g{t)y and X = ip, where / and g run through the set of smooth functions of t. The physically 
most important of these characteristics are A = 1, A = y and \ = ip, which are associated with 
the conservation of circulation, x- momentum and energy. Any zero-order characteristic of Eq. (1) 
is a characteristic of Eq. (18). Indeed, denoting L := C,t + ipxCy ~ i^yCx + l^ipx — D we derive that 

fL = (^fip^t + my + fPi^ - i^f^x^^ + By (^fi^yt - mx - i^fDy^^ , 

gyL = (^gyipxt + gyipCy - ^{i^yf + Qvl^i^ - T^gy^x^^-^ 

( AC^ AC^\ 

+ Dy ( ^yV's/t - gi'y - gyPCx + gipipxy - vgyDy—^ + vg—^ j , 

iPL = Bt (^-i(V^)2^ + D,. (^V'V'x.t + \^\y + ^V'' - l^i^^x^ + lyi'x^ - l^B^C'^ 

f 1 2 AC^ vA 
+ Dy ( V'^s/t - Cx - i"iljT>y—^ + yijjy— z^DyC 1 • 

If we grant that the vorticity equation coupled with some diffusive term possesses a smaller 
number of conservation laws (e.g. owing to the special physical properties of this diffusion), we 
can use a simpler form for the expression D. For example, the differential function D = v/S.C,^ 
leads to a parameterization which is invariant with respect to the entire pseudogroup Gi and 
possesses conservation laws with the characteristics A = /(t), A = g{t)y for arbitrary values of 
the smooth parameter-functions / and g. 

The parameterization (18) demonstrates the feasibility of combining invariant and conserva- 
tive properties of closure schemes. This possibility is important for two obvious reasons. Firstly, 
conservation laws incorporate relevant physical information that is worth being preserved by 
a parameterization scheme. Secondly, from the point of view of constructing parameterization 
schemes, the requirement of preserving both symmetries and conservation laws leads to a more 
specific class of schemes than considering either only symmetries or only conservation laws. The 
additional narrowing of the class of admitted schemes using geometric constraints can then help 
to reduce the number of schemes that must be tested numerically so as to find the optimal 
parameterization for a given process. 



9 Conclusion and discussion 

The differential invariants of the Lie symmetry pseudogroup Gi of the barotropic vorticity equa- 
tion on the beta-plane are computed using the technique of moving frames for Lie pseudogroups. 
A basis of these differential invariants along with the associated operators of invariant differ- 
entiation is established. Together, they serve to completely describe the algebra of differential 
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invariants of Gi. Although differential invariants have many applications (such as the inte- 
gration of ordinary differential equations [35], computation of so-called differentially invariant 
solutions [22, 39], the construction of invariant numerical discretization schemes [17], etc.), in 
the paper we focus on their usage in the construction of invariant closure schemes or, perhaps 
more generally, invariant diffusion terms for the averaged vorticity equation. This is one of 
the two general methods proposed in [41] to derive parameterization schemes with symmetry 
properties. As an alternative to the direct usage of elementary differential invariants that can be 
build together to yield invariant closure schemes, we propose the method of invariantization of 
existing parameterization schemes. This method is along the line of invariantization of existing 
discretization schemes as introduced in [25, 26]. Though straightforward to apply, a potential 
disadvantage of this method is that the result depends on the particular choice of the moving 
frame and therefore does not lead to a unique invariant counterpart of an existing noninvariant 
scheme. As a consequence, it might be necessary to modify invariantized closure schemes and 
to test different invariantizations in order to devise physically valuable closures. 

The differential invariants derived are used to construct invariant hyper diffusion terms in 
order to model the behavior of two-dimensional freely decaying turbulence. The resulting en- 
ergy spectra exhibit a region of approximate slope which is the theoretically derived shape 
for the postulated enstrophy inertial range. It should be stressed, though, that the obtained 
energy spectrum should be taken with a pinch of salt. Since the derivation of the theoretical 
form of the spectra in [4, 28] it has been tried in numerous studies to obtain these spectra 
in numerical simulations. Although results often vary, spectra are found with a steeper slope 
than the predicted curve as described in [6, 7, 13, 29, 31, 48]. It seems to be generally 
agreed today that the presence of the stable coherent vortices, which is the main feature of 
two-dimensional turbulence, has a strong impact on the derived energy spectra. This holds 
both in the case of turbulence on the f-plane and on the beta-plane. The introduction of an 
invariant hyperdiffusion-like term certainly complicates the situation as diffusion then is coupled 
nonlinearly to the vorticity equation. On the other hand, it was indicated that the presence 
of the beta-term in the vorticity equation allows for a nonlocal transfer of anisotropy from the 
larger to the smaller scales [31]. A nonlinear diffusion term has the potential to support such a 
nonlocal scale interaction and thereby serves as a potential parameterization scheme for numer- 
ical models. It should be stressed in this context that in all the simulations we have carried out, 
the rate of energy dissipation was significantly lower than using classical hyperdiffusion even in 
quite low-resolution numerical experiments. 

Apart from the discussion above, the possibility of constructing hyperdiffusion-like enstrophy 
sink terms that lead to scale invariant energy spectra seems to be a valuable property for 
itself. It is precisely the scale invariance of the Euler equations that is used to predict the 
behavior of two-dimensional turbulence in the inertial range and therefore the availability of 
dissipative versions of the vorticity equation having the same invariance properties as the inviscid 
vorticity equation might be a general advantage. Heuristically, one can expect that an invariant 
closure scheme should be better adapted for the problem of reproducing features that have 
been derived using symmetries (as the isotropic energy spectrum), similarly as an invariant 
discretization scheme often reproduces better invariant exact solutions of a differential equation 
than noninvariant discretization schemes [44]. This assumption is supported by the proved 
relevance of Lie symmetries in turbulence theory [34] . The results obtained in the present paper 
do not contradict this assumption, keeping in mind especially that the premises invoked to obtain 
the theoretical form of the spectra are at present under revision. In this context, it should again 
be stressed that there is a multitude of invariant parameterization schemes or invariant diffusion 
terms that can be coupled to the vorticity equation on the beta-plane. The fact that already the 
simplest invariantized version (16) of the hyperdiffusion term (which has obvious weaknesses) 
shows quite good properties in the course of our numerical tests is a motivating result which 
is worth pointing out. Nevertheless, in order to verify and better assess the ability of invariant 
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hyperdiffusion schemes to model turbulence on the beta-plane, further theoretical and numerical 
studies must be carried out. 

Note that the possibility of constructing diffusion schemes that are invariant under the en- 
tire maximal Lie invariance group is a particular feature of the beta-plane vorticity equation, 
which is most probably not realizable as straightforward for vorticity dynamics on the f-plane. 
The complication with the latter model is that the corresponding maximal Lie invariance pseu- 
dogroup Go is even wider than Gi. This makes it much harder to derive reasonably simple 
closure schemes that are invariant under the entire pseudogroup Go, see the discussion in [41], 
where a generating set of differential invariants of Go and a complete set of its independent 
operators of invariant differentiation are determined. A possible remedy for this complication 
is to consider closure schemes that are invariant only under certain subgroups of the maximal 
Lie invariance pseudogroup of the f-plane equation. As highlighted in the present paper, the 
selection of such subgroups can be justified for physical reasons when boundaries come into play. 

Another novel feature of the present paper is the explicit inclusion of conservation laws in 
invariant closure schemes. The chance of constructing such conservative invariant parameter- 
ization schemes is of obvious physical relevance. For physical processes that do not violate 
particular conservation laws, it is natural to require the associated parameterization to be also 
conservative. It was demonstrated in the paper for the vorticity equation on the beta-plane 
that the concepts of invariant and conservative parameterization schemes can be united to yield 
closure ansatzes that preserve both all the symmetries and certain conservation laws of this 
equation. The construction of further invariant conservative closure schemes as well as their 
exhaustive testing will be a next major challenge in the application of ideas of group analysis to 
the parameterization problem. 
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